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ABSTRACT 


The  purpose  of  this  thesis  is  to  examine  the  feasibility  of  using  an  emerging 
technique,  called  transformation  optics  (TO),  for  designing  materials  to  be  used 
as  a  defense  against  directed  energy  weapons  for  satellites.  In  order  to  do  this,  a 
method  of  determining  the  effectiveness  of  TO  against  high-intensity  fields  must 
be  demonstrated.  These  high-intensity  fields  will  cause  a  nonlinear  response  in 
the  material  and  it  is  this  nonlinear  response  that  will  be  studied.  TO  has  been 
shown  to  be  effective  when  dealing  with  lower  intensity  fields  and  thus  linear 
responses  in  matter.  [1]  This  thesis  will  attempt  to  model  the  nonlinear  response 
and  solve  for  the  fields  due  to  this  response. 

The  fields  induced  by  the  nonlinear  response  are  considered  an  error  field. 
To  solve  for  the  error  field,  a  method  to  model  the  nonlinear  response  will  be 
derived  using  Miller’s  Rule.  Stemming  from  the  Lorentz-Drude  model  of 
polarization,  Miller’s  Rule  serves  as  a  model  of  the  nonlinear  response  and  has 
been  shown  experimentally  to  be  approximately  true.  [2]  Once  the  nonlinear 
response  has  been  found,  the  error  can  be  analyzed  as  an  electrostatic  problem 
to  determine  if  the  polarization  or  magnetization  induces  a  field  within  the  cloaked 
area. 
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I.  INTRODUCTION 


A.  THE  PROBLEM  DEFINED 

The  United  States  military’s  reliance  on  satellites  is  well  known.  Whether  it 
is  intelligence  collection,  communication,  or  precision  navigation  and  timing,  the 
United  States  relies  heavily  on  its  use  of  satellites  in  space  to  conduct  effective 
and  efficient  operations  around  the  globe.  As  potential  adversaries  look  to  exploit 
this  reliance,  antisatellite  weapons  in  the  form  of  directed  energy  beams  are  a 
potential  future  threat  to  our  satellites  that  cannot  be  overlooked. 

When  directed  energy  weapons  are  mentioned,  often  the  first  things  that 
come  to  mind  are  the  laser  beams  of  science  fiction.  In  these  stories,  the  usual 
defense  is  some  sort  of  “energy  shield.”  While  it  may  be  prominent  only  in 
science  fiction,  shielding  from  electromagnetic  waves  is  possible  and  has  been 
demonstrated  using  cloaking  techniques.  [1]  Cloaking,  in  this  sense,  is  the  ability 
of  certain  designed  materials  to  redirect  electromagnetic  waves  around  a  region 
in  space,  thus  cloaking  the  object.  While  the  end  goal  of  cloaking  is  invisibility,  for 
a  directed  energy  weapon  simply  redirecting  the  wave  is  sufficient.  The  fields 
(and  thus  energy)  involved  with  current  demonstrations  so  far  has  been  weak, 
dealing  strictly  in  the  linear  realm.  [3]  With  directed  energy  weapons,  it  will  be 
necessary  to  look  at  nonlinear  effects.  This  thesis  will  deal  with  modeling 
nonlinear  responses  of  materials  specifically  designed  to  shield  against  directed 
energy  weapons. 

Within  the  last  decade,  a  technique  for  designing  materials  that  exhibit  the 
desired  cloaking  response  has  been  developed  called  transformation  optics.  The 
field  of  transformation  optics  has  introduced  new  methods  of  analysis  and  design 
of  electromagnetic  materials.  The  method  simplifies  many  problems  in 
electrodynamics  by  converting  what  would  be  a  complicated  analytical  solution 
into  a  geometric  coordinate  transformation.  [4]  The  trade-off  is  that  the  material 
needed  to  implement  this  transformation  in  real  space  is  often  complicated  to 


1 


manufacture.  This  hurdle  has  been  overcome  somewhat  in  the  last  decade  and 
transformation  optics  is  becoming  a  viable  design  technique  for  a  whole  host  of 
new  materials.  [1] 

Typically,  transformation  optic  solutions  are  implemented  with  the  use  of 
metamaterials.  Metamaterials  are  “artificially  structured  media  whose  effective 
constitutive  parameters  can  be  adjusted  over  a  wide  range  of  parameters 
through  judicious  design.”  [5]  For  the  purpose  of  designing  a  shield  for  a 
satellite,  weight  is  one  obvious  concern.  Luckily  these  metamaterials  can  be 
manufactured  from  lightweight  alloys  and  do  not  require  large  scale  structure.  [1] 
In  fact,  the  structures  are  often  sub-wavelength  in  size.  In  this  thesis,  the 
effectiveness  of  the  technique  of  transformation  optics  for  designing 
metamaterials  for  use  in  space  is  evaluated  and  if  the  nonlinear  effects  of  the 
high  intensity  fields  associated  with  directed  energy  weapons  can  be 
approximated  through  Miller’s  rule. 

Miller’s  rule  is  a  technique  to  predict  the  second  order  nonlinear 
susceptibility  of  a  material  in  a  consistent  way,  in  close  analogy  to  a  power  or 
Taylor  series  expansion.  Based  on  the  first  order  linear  response,  an 
approximation  is  made  for  the  second  order  response.  Once  the  second  order 
susceptibility  is  known,  the  second  order  polarization  is  calculated.  The  nonlinear 
response  is  then  analyzed  as  an  electrostatic  problem  to  determine  if  the 
polarization  or  magnetization  induces  a  field  within  the  cloaked  area. 

B.  ANTISATELLITE  (ASAT)  TECHNOLOGY 

1.  Definition 

Antisatellite  (ASAT)  technology  has  gained  attention  in  recent  years  due  to 
two  high  profile  tests  of  ASAT  weapon  systems  by  the  Chinese  and  the  United 
States.  The  sanctuary  provided  by  the  altitude  and  speed  of  low  earth  orbit  (LEO) 
has  been  lost.  The  United  States’  reliance  on  satellite  systems  cannot  be  over¬ 
stated,  and  the  subsequent  threat  posed  by  antisatellite  weapon  systems  cannot 
be  ignored. 
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The  purpose  of  a  weapon  system  is  to  deliver  a  destructive  amount  of 
energy  to  a  target  to  damage  or  destroy  it.  Conventionally,  a  projectile  fired  at 
great  velocity,  an  explosive  charge,  or  a  combination  of  both  transfers  the  energy 
to  the  target.  The  larger  the  energy  transfer,  the  greater  the  destruction.  Directed 
energy  weapon  systems  (DEWS)  use  electromagnetic  waves  to  transfer  energy. 
They  focus  or  direct  large  amounts  of  energy  onto  a  relatively  small  spot.  By 
concentrating  the  energy,  DEWS  require  less  total  energy  to  inflict  damage. 

a.  Conventional  AS  AT 

The  two  ASAT  systems  recently  tested  are  kinetic  in  nature, 
requiring  launch  by  the  host  country  and  intercept  capability  by  the  ASAT.  This 
requirement  of  highly  advanced  technology  to  build  an  ASAT  has  served  to  limit 
the  proliferation  of  antisatellite  systems.  This  could  change  in  the  future.  “Robert 
Joseph,  the  State  Department’s  point  man  for  arms  control  and  international 
security,  said  other  nations  and  possibly  terrorist  groups  were  ‘acquiring 
capabilities  to  counter,  attack  and  defeat  U.S.  space  systems.’  “  [6] 

In  addition  to  the  technical  hurdles  to  overcome,  the  destruction  of 
a  satellite  in  LEO  poses  risks  to  other  satellites,  including  manned  space  flight, 
within  the  LEO  belt.  International  ramifications  from  the  Chinese  test  serve  to 
support  this.  “A  variety  of  regional  and  extra-regional  states  have  expressed 
concern  about  a  suspected,  although  not  publicly  confirmed,  antisatellite  weapon 
(ASAT)  test  by  China  on  1 1  January.”  [7]  The  debris  cloud  formed  from  the 
satellite  and  exploding  missile  is  a  problem  for  all  nations  wishing  to  utilize  the 
LEO  belt.  “According  to  David  Wright  of  the  Cambridge,  Mass. -based  Union  of 
Concerned  Scientists,  the  satellite  pulverized  by  China  could  have  broken  into 
nearly  40,000  fragments  from  1  to  10  centimeters  (a  half-inch  to  4  inches)  in  size, 
roughly  half  of  which  would  stay  in  orbit  for  more  than  a  decade.”  [8]  If  an  ASAT 
system  could  be  developed  that  limits  debris  and  requires  simpler  technology  to 
operate,  then  the  above  obstacles  could  be  overcome. 
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b. 


Directed  Energy 


One  class  of  such  weapons  are  directed  energy  weapons.  The  best 
known  example  of  a  directed  energy  weapon  would  be  a  laser,  but  other  forms 
exist,  including  directed  microwave  weapons.  For  the  reasons  mentioned  above, 
the  U.S.  must  develop  a  robust  counter-directed  energy  weapon  (CDEW) 
capability  against  directed  energy  weapons. 

Directed  energy  weapons  offer  two  distinct  advantages  over 
conventional,  kinetic  antisatellite  weapons.  The  first  is  the  ease  of  delivering 
energy  to  LEO.  For  a  conventional  system,  the  host  nation  needs  to  have 
developed  a  launch  and  rendezvous  capability  of  extremely  high  precision.  This 
is  currently  beyond  the  reach  of  many  nations  and  not  something  easily  acquired. 
The  second  advantage  is  the  destruction  of  the  satellite,  or  at  least  its  capability, 
without  the  debris  cloud,  and  thus  without  the  international  ramifications.  A 
directed  energy  weapon  could  simply  cause  enough  damage  by  burning  through 
or  creating  fields  within  the  electronics  to  damage  them.  The  satellite  would 
become  inoperable  yet  remain  in  the  same  orbit.  Directed  energy  weapons  are 
currently  being  developed  for  these  reasons.  “The  latest  report  follow  claims  in 
September,  reported  by  Defense  News,  that  China  was  aiming  high-powered, 
ground-based  lasers  at  U.S.  spy  satellites  —  apparently  to  test  whether  sensors 
on  the  satellites  could  be  blinded.”  [6] 

One  defensive  technology  that  could  serve  to  protect  U.S.  satellites 
by  bending  electromagnetic  fields  away  from  a  region  of  space  is  TO.  The 
bending  or  cloaking  effect  is  typically  implemented  through  metamaterials.  This 
cloaking  would  render  the  satellite  invisible  to  the  directed  energy,  causing  it  to 
pass  harmlessly  around  or  away  from  the  asset  requiring  protection.  While  still 
highly  theoretical,  examples  of  cloaking  a  cylinder  in  two  dimensions  from 
microwave  radiation  have  been  demonstrated.  [1]  The  fields  involved  in  this 
experiment  were  of  low  intensity  however,  and  did  not  demonstrate  an  effect  for 
high  intensity  fields.  For  cloaking  to  be  effective  at  high  intensities,  nonlinear 
effects  need  to  be  taken  into  account. 
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II.  ELECTRODYNAMICS  OF  MATTER 


A.  MAXWELL’S  EQUATIONS  IN  MATTER 

1.  Introduction 

In  order  to  understand  transformation  optics,  it  is  necessary  to  understand 
the  behavior  of  electromagnetic  waves  in  media.  The  reader  is  assumed  to  be 
familiar  with  electrostatics,  Maxwell’s  equations,  and  general  electrodynamics. 
This  section  is  included  for  completeness  but  does  not  serve  as  a  complete 
derivation  or  proof,  rather  as  a  review  and  reference  for  subsequent  sections. 
Brau’s  Modern  Problems  in  Electrodynamics  and  Mill’s  Nonlinear  Optics  are  both 
excellent  sources  for  further  study  of  the  material  in  this  section.  [2]  [9] 

2.  Maxwell’s  Equations  in  Matter 

Maxwell’s  equations  are 


V  -E  =  — 

e„ 

(2.1) 

V  B  =  0 

(2.2) 

VxE  =  -  — 
dt 

(2.3) 

dE 

M,,\  +  A, A  — 
dt 

(2.4) 

Here  E  is  the  electric  field,  B  the  magnetic  induction  field,  p  the  charge  density, 
and  j  the  current  density.  Maxwell’s  equations  are  true  everywhere,  to  include 
inside  media.  They  are,  however,  difficult  to  apply  due  to  a  lack  of  knowledge  of 
the  microscopic  structure  of  most  media.  Nonetheless,  field  averages  can  be 
examined. 

When  a  material  is  immersed  in  an  electric  field,  it  becomes  polarized.  An 
additional  field  P,  called  the  “polarization  field,”  is  then  induced  due  to  the 

5 


molecular  interactions  with  E.  It  can  be  shown  that  P  is  the  dipole  moment  per 
unit  volume.  In  this  case,  (2.1 )  becomes 


V  -(E  +— P)  =  — 

(2.5) 

Define  a  new  field  D,  called  the  “displacement  field” 

D  s  £oE  +  P 

(2.6) 

In  this  case  (2.1)  becomes 

V  D  =  p 

(2.7) 

In  a  similar  way,  matter  responds  to  an  imposed  magnetic  field  through 

the  magnetization  response  M. 

Define  the  magnetic  field  H  as 

H  =  — B-M 

Mo 

(2.7) 

in  which  case,  (2.4)  becomes 

VxH  =  j  +  — 
dt 

(2.8) 

Combining  the  above  results,  Maxwell’s  equations  in  matter  are 

V-D  =  p 

(2.9) 

V  -B  =  0 

(2.10) 

„  „  SB 

VxE  = - 

dt 

(2.11) 

5D 

VxH  =  j  +  — 
dt 

(2.12) 

Of  primary  concern  is  the  polarization  field  P  induced  by  high  intensity  electric 
fields.  The  next  section  examines  P  more  closely,  in  particular  its  response  to  an 
applied  high  intensity  electric  field. 
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B.  POLARIZATION  RESPONSE  TO  LARGE  ELECTRIC  FIELDS 


1.  Derivation  [9] 

If  the  applied  electric  field  is  of  sufficient  strength  to  invoke  a  measurable 
nonlinear  response  in  the  material  yet  not  strong  enough  as  to  overcome  the 
coulomb  attraction  between  the  electron  and  nucleus,  then  the  polarization  can 
be  expanded  as  a  Taylor  series  in  powers  of  the  macroscopic  electric  field  such 
that 


^(r.O^+Z 


f  dP  A 


8Ea 


I Eg+- 
2! 


-Z 

1 1  i—t 


Pr 


dEpdEy  j 


\EPE7+.., 


n: 


!  (n-r)l 


(2.13) 


Here  Pa(r,r)is  the  ath  Cartesian  component  of  the  dipole  moment  per  unit 
volume  with  a  ranging  over  x,  y,  and  z. 

In  most  cases,  the  first  term  of  equation  (2.14)  vanishes  identically.  If  it  did 
not,  the  material  would  have  an  intrinsic  polarization,  analogous  to  a 
ferromagnetic  material.  This  is  equivalent  to  stating  that  any  polarization  present 
is  due  to  the  applied  electric  field.  With  this  being  the  case,  equation  (2.14) 
becomes 


(  flp  |  1 

PA r,t)  =  y  — ^  EJr,t)+—Y 
“  ^  dEn  \  p  2\j? 


P  J 


s2p 

dEpdErj 


\EJr,t)E(r,t)  +  ... 


(2.14) 


Based  on  (2.15),  the  susceptibilities  of  the  material  are  defined  as  the  partial 
derivative  terms.  In  other  words, 


y(1)  = 


dE„ 


apy  dEpdEy 


(2.15) 

(2.16) 


are  the  first  order  and  second  order  susceptibilities  respectively.  The  quasi-static 
limit  will  be  used  in  this  analysis.  It  will  be  assumed  that  a  steady  state  condition 
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has  been  achieved  to  remove  the  time  dependence  from  (2.15).  Given  the 
above,  the  polarization  is  now 

P.  O')  =  VO  +  Z*  (r)  + ...  (2.17) 

As  can  be  seen,  transformation  optics  requires  tensor  analysis.  The 

metamaterials  themselves  are  intrinsically  inhomogeneous  and  anisotropic, 

which  requires  that  the  susceptibilities  be  tensors. 

The  polarization  is  now  separated  into  a  linear  and  nonlinear  contribution. 

Pa{r)  =  P^L\r)  +  Pr\r)  (2-18) 

From  (2.18),  the  linear  response  is 

p?'  <2'19> 

and  the  nonlinear  response  is 

p-'"’(r)  =  Sg,2'S^(r)£,(r)  +  ...|  (2.20) 

For  this  analysis,  the  second  order  response  will  be  examined  with  higher  order 
contributions  ignored.  E.g.  (2.21)  describes  the  nonlinear  polarization.  The  next 
step  is  to  determine  the  nonlinear  susceptibility  using  Miller’s  rule. 

C.  LORENTZ-DRUDE  MODEL  OF  POLARIZATION  AND  MILLER’S  RULE 
DERIVATION  FOR  SECOND  ORDER  SUSCEPTIBILITY  [2] 

Miller’s  rule  is  a  method  that  approximates  the  second  order  susceptibility 
in  terms  of  the  first  order  susceptibility.  A  simple  and  elegant  derivation  of  Miller’s 
Rule,  paraphrased  in  this  section  from  Brau’s  textbook  Modern  Problems  in 
Electrodynamics,  follows  from  the  Lorentz-Drude  Model  for  the  polarization  of  the 
atom[2].  In  this  model,  the  electron  is  harmonically  bound  to  the  nucleus.  In  many 
cases,  it  may  be  difficult  to  derive  the  first  order  susceptibility  of  a  material,  but 
with  transformation  optics,  the  derivation  is  simple.  As  will  be  shown  later,  the 
first  order  susceptibility  is  simply  a  function  of  the  transformation  matrix  used  in 
transformation  optics.  The  derivation  of  the  second  order  susceptibility  begins 
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with  a  simple  linear  isotropic  case,  extends  to  the  nonlinear  case,  and  finally  to 
anisotropic  materials. 


1.  ISOTROPIC  MATERIALS 
a.  Linear  Materials 

In  the  presence  of  an  electric  field  E,  the  displacement  r  of  the  atom 
in  the  nonrelativistic  limit  is 


dr  2  d  t~i 

+  y  —  +  co~r  =—  E 
dt  '  m 


(2.21) 


Here  q  is  the  charge,  m  the  mass,  y  the  damping  coefficient,  and  cor  the  resonant 
frequency  of  the  system.  The  polarization  per  unit  volume  is  the  dipole  moment 
qr  times  the  electron  density  ne,  or  P  =  negr.  This  leads  to  the  equation 


where 


d  P  dP  i  2 

-zr+y—+a)-rP  =  £oa)  pE 


dt 1 


dt 


(2.22) 


0)„ 


(2.23) 


is  the  plasma  frequency  of  the  material.  Taking  the  Fourier  transform  of  both 
sides  of  this  equation,  keeping  in  mind  that  d/dt  transforms  to  — /co, 

( c o 2  -co2  - iyco)P  =  s0 co2E  (2.24) 

Here  the  tilde  indicates  the  transformed,  or  frequency  domain,  fields. 

Comparing  the  above  equation  with  the  equation  for  polarization 

P  =  eoZE  (2.25) 


it  can  be  seen  that  the  dielectric  susceptibility  is 
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z(®)  = 


(2.26) 


ojr  -  or  -  iyo) 

b.  Nonlinear  Materials  [2] 

Taking  this  result  and  extending  it  to  the  nonlinear  realm  requires 
an  additional  term  in  the  inharmonic  oscillator  equation.  The  new  equation 
describing  the  electron  bound  to  the  nucleus  is 

^  +  y  —  +  co2l.P  +  K(2)P2  =tjE  (2.27) 

dt2  dt 


where  K(2)  represents  the  nonlinear  contribution  to  the  restoring  force.  The 
parameter  ii~\  is  introduced  to  act  as  an  ordering  parameter. 

Assume  that  r\  is  sufficiently  small  such  that  the  solution  can  be 
written  as  a  power  series  having  the  form 

P(t)  =  riP(l)(t)  +  77  2P(2\t)  + ...  (2.28) 

Take  the  Fourier  transform  of  this  equation  to  obtain 

P(oj)  =  T]Pm(co)  +  77  2Pr-\co)  +  ...  (2.29) 


In  order  to  compute  the  polarization  now,  plug  the  above  equation  into  (2.28)  and 
collect  like  terms  of  r\.  To  first  order  in  r\ 


d2Pm  dP(l) 

- 7 — ^  Y  — 

dt'  dt 


+  (o2  rP 


(1)  =E 


(2.30) 


The  Fourier  transform  is  again  taken  and  compared  to  the  answer 
of  the  equation  for  polarization  (2.26). 
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lP"\a o )  =  -g-  ,  g:  .  E(m) 
s0  cor  -co  -  iyco 


(2.31) 


From  this,  it  is  seen  that  the  first  order,  linear  susceptibility  is 


s0  cor  -  co  -  iyco 


(2.32) 


Taking  the  exact  same  steps,  but  this  time  to  the  second  order  of  r\ 
the  equation  obtained  is 


d2P  dP [ 

^4-  +  y—  +  co2rP(2)+K(2\Pm)2  =  0 

dt2  dt 


(2.33) 


Again  taking  the  Fourier  transform  and  using  (2.32),  it  is  found  that 


(co2  -co2  -  iyco) P(2\co)  =  - 


co 

==  dt  eUot(Pa\t)f 

L 


Its- (2)  co  co  co 

-  0  3/2  2  \dt\ dco'f 

V  ^  )  V  —oo  —00  —00 


(2.34) 


The  integral  over  t  is  a  delta  function.  Using  (2.33),  the  nonlinear  polarization  is 


rj2P(2)(co)  =  £0\  dco'\ dco"x(2\co\co")E(co')E(co")S(co-co'-co")  (2.35) 


where  the  second  order  susceptibility  is 
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/(2)(©^©'0  =  -^^/)(©'+©")/)(®')/(1)(®") 


(2.36) 


This  is  Miller’s  rule. 


2.  ANISOTROPIC  MATERIALS 
a.  Nonlinear  Materials  [2] 

The  above  derivation  for  Miller’s  rule  applies  to  isotropic  materials. 
Metamaterials  are  by  their  nature  anisotropic.  In  order  to  derive  the  anisotropic 
Miller’s  Rule,  proceed  as  before  but  now  expand  into  three  dimensions  by 
including  the  appropriate  indices.  The  following  derivation  includes  a 
nonstandard  notation,  utilizing  a  semicolon  to  denote  an  index  which  does  not 
exhibit  symmetry.  The  indices  on  the  right  of  the  semicolon  are  permutable;  the 
index  on  the  left  is  not.  This  semicolon  can  be  ignored  for  the  most  part,  but 
offers  a  way  to  simplify  the  calculations.  In  other  words  %.  k  =  xrk, . 

Assume  the  solution  has  the  form 


pI{t)  =  nP^\t)  +  v2^2\t)  +  ... 


(2.37) 


with  a  Fourier  transform 


P(co)  =  i]Pf1)(co)  +  rfPf2)(co)  +  ... 


(2.38) 


Inserting  equation  (2.38)  into  the  equation  of  motion  to  first  order  r|  obtains 


2nd) 


dlp: 


dt2 


V  dp< 


(1) 


+ 


z 


CO 


\JpJ 


(!)  _  £(1) 


(2.39) 
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The  Fourier  transform  of  (2.40)  is 


where 


I 


(1)  _ 


=E 


K.  -CO2 . 

i;j  >,j 


-iooy,j 


■  co2d:: 


The  solution  to  this  equation  is 


^<1)=foX^ 


where  the  linear  susceptibility  tensor 


X(X)=—k:1 

kj 

bn 


is  also  the  inverse  of  the  tensor  Ar...such  that 


£z$Kj*  =  >/4 

j 


Using  the  same  arguments  as  the  previous  cases,  the 
term  in  the  polarization  is 


(2.40) 

(2.41) 

(2.42) 

(2.43) 

(2.44) 

quadratic 
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00  00 

T12Pf2)((o)  =  £0YJ  J  deaf  J  dofX^k  (©',  of)  Ej  (co')  Ek  (g)')S(g)  -co'-  of)  (2.45) 

j)k  —00  —00 


where  the  second  order  susceptibility  tensor  is 


y(2) 
A-  i ;  jk 


s~K[ 

_  V  0  pw 

\G)  ,  0)  )  2_j  r~ — 

p,q,r 


xZ  +  a")Zaj  (a')zZ  (®") 


CD  i 


(2.46) 


Equation  (2.47)  is  the  Miller’s  rule  determined  second  order 
nonlinear  susceptibility.  Since  the  quasi-static  case  is  being  examined,  all 
frequency  dependence  is  removed. 


(2.47) 


Additionally,  it  can  be  shown  [2]  that  there  exists  symmetry  across  the  last  two 
indices  of  the  second  order  susceptibility,  but  not  the  first.  This  is  a  cross  check 
of  the  calculation. 


While  the  first  order  susceptibility  is  a  unitless  value,  the  second 
order  susceptibility  has  the  units  of  meters  per  Volt,  or  m/V,  in  the  SI  system. 
Also  of  note  is  the  second  order  susceptibility  involves  the  product  of  three  first 
order  susceptibilities.  It  then  makes  sense  that  a  material  with  a  strong  first  order 
susceptibility  will  have  a  relatively  strong  second  order  susceptibility  as  well.  [2] 


b.  Miller’s  coefficient 

As  seen  in  (2.48),  there  is  a  collection  of  terms  that  must  be 
accounted  for  when  calculating  the  nonlinear  susceptibility 

£2K{2) 

(2.47) 

v2tzt7 
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This  collection  of  terms,  sometimes  referred  to  as  Miller’s  coefficient  when  taken 
together,  has  been  found  to  be  approximately  equal  to  4.52  x  10-2  m2C'1  for  a 
wide  range  of  materials  studied.  [10]  The  reason  for  this  is  not  well  understood, 
but  it  is  found  to  be  experimentally  true  to  within  an  order  of  magnitude.  [10]  For 
anisotropic  materials,  this  coefficient  becomes  a  tensor,  further  complicating  its 
calculation.  Referring  back  to  (2.48)  and  removing  the  directional  dependence  for 
simplicity,  the  coefficient  is  defined  to  be 


4k{2) 


(2.47) 


Here,  77  is  the  ordering  parameter,  but  it  can  also  be  shown  to  equal 


m 


(2.47) 


where  N  is  the  number  of  oscillators,  q  the  charge,  and  m  the  mass.  One  can 
see  that  to  calculate  this  value  directly  for  an  anisotropic  metamaterial  would  be 
extremely  difficult.  It  will  therefore  be  set  to  its  empirical  value. 
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III.  TRANSFORMATION  OPTICS 


A.  BACKGROUND 

TO  uses  coordinate  transformations  to  determine  material  specifications 
that  control  electromagnetic  fields  in  interesting  and  useful  ways.  [11]  By 
applying  a  coordinate  transformation  to  the  constitutive  relations,  a  material  can 
be  designed  that  mimics  free  space  yet  possesses  a  curved  geometry.  It  is 
important  to  keep  in  mind  that  in  the  general  case  the  constitutive  relations  are 
tensor  quantities. 

TO  allows  the  derivation  of  a  linear  constitutive  relation  directly  from 
a  desired  trajectory  of  light  -  in  a  completely  algebraic  way.  The 
key  ingredient  of  TO  is  to  define  the  trajectory  in  the  specified 
medium  as  a  result  of  a  transformation  applied  on  the  trajectory  in 
free  space.  Unlike  normal  coordinate  transformations,  under  which 
physics  is  invariant,  the  field  lines  are  considered  to  be  attached  to 
the  coordinates  during  the  transformation,  thereby  changing  the 
physics  and  deriving  the  constitutive  law  of  the  desired  medium. 

[12] 

Analogous  to  general  relativity,  the  medium  in  effect  fools  the  electromagnetic 
wave  into  believing  it  is  propagating  in  free-space. 

B.  LINEAR  FORMULATION 

Smith,  Kundtz,  and  Pendry  give  a  succinct  explanation  of  the  concepts 
behind  transformation  optics. 

It  has  been  recognized  for  some  time  that  Maxwell’s  equations  can 
be  written  in  a  form-invariant  manner  under  coordinate 
transformations,  such  that  only  the  permittivity  and  permeability 
tensors  are  modified.  With  the  coordinate  transformation  applied  to 
the  constitutive  parameters,  electromagnetic  waves  in  one 
coordinate  system  can  be  described  as  if  propagating  in  a  different 
coordinate  system.  [5] 

Therefore,  the  first  task  will  be  to  write  Maxwell’s  equations  in  this  form-invariant 
manner. 
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1.  Arbitrary  Coordinates  [4] 

It  is  first  necessary  to  introduce  tensor  notation  and  its  associated 
concepts.  This  section  is  not  intended  to  be  a  full  discussion  on  the  topic,  but 
rather  it  serves  as  a  review  of  the  concepts  necessary  for  TO.  For  a  more 
thorough  handling  of  the  differential  geometry  involved,  see  Leonhardt  and 
Philbin’s  article  “Transformation  Optics  and  the  Geometry  of  Light”  included  in 
the  bibliography.  The  discussion  below  borrows  heavily  from  section  three  of  the 
article.  [4] 


a.  Einstein  Summation  Convention 

The  Einstein  summation  convention  is  defined  as  the  summation 
over  any  repeated  index,  whether  it  be  in  a  subscript  or  superscript,  over  its 
entire  range.  As  an  example 

aixi  =  OjXj  +  a2x2  +  a3x3  (3 . 1 ) 

Additionally,  the  location  of  the  indices  will  have  meaning.  A  superscript  will  be 
taken  to  be  a  contravariant  vector  (column  vector)  while  a  subscript  is  a  covariant 
vector  (row  vector). 

b.  Coordinate  Transformations 

TO  deals  with  general,  curvilinear  coordinates  and  it  is  necessary  to 
have  a  way  of  expressing  a  generalized  coordinate  system  and  performing  an 
arbitrary  transformation  to  another  set  of  coordinates.  Define  one  coordinate 

system  |x',/  =  l,2,3|  that  is  transformed  to  another  system,  differentiated  with  a 

prime jx'',/’  =  1,2,3] .  Here  x’"  =  [x']'  or  the  i-th  component  in  the  primed 

coordinate  system.  The  differentials  of  this  coordinate  system  are  related  by  the 
chain  rule. 

dx'  =  —v  dx’  ,  dx’  =  —  dxl  (3 . 1 ) 

dx‘  dx 
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In  order  to  carry  out  the  transformation,  define  a  standard 
transformation  matrix 


II 

■<f 

(3.1) 

such  that 

dx'  =  A  lvdxl  =  A'., A  ljdx] 

(3.2) 

and 

dx'  =A  '.dx1  =  A'.A '  vdxJ 

i  1  J 

(3.3) 

Taken  together,  the  above  two  equations  imply 

KK  =  S) 

AfA  'r  =  S'r 

(3.4) 

Where  <5'  is  the  Kronecker  delta.  In  other  words,  A‘,  and  A''are  inverses  of  each 

J  i  i 

other. 


c.  The  Metric  Tensor 

Using  the  above  notation,  any  arbitrary  coordinate  system  can  be 
defined.  One  thing  that  must  remain  invariant  however  is  the  distance  between 
two  points,  despite  the  coordinate  system  used  to  represent  those  points.  In 
mathematics,  the  Pythagorean  theorem  is  a  relation  in  Euclidean  geometry 

ds2  =  dx2  +  dy2  +  dz2  (3-5) 

which  is  expressed  in  the  new  notation  as 

ds2  =  Sjdx'dx1  (3.6) 

where  AJs  again  the  Kronecker  delta.  For  arbitrary  coordinates  however,  the 
distance  between  two  points  is 

ds1  =  gijdxldxi  (3.7) 
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The  new  term  in  the  above  equation,  g  is  known  as  the  metric 
tensor.  It  is  a  symmetric  tensor(g;>.  =£,,),  which  allows  the  calculation  of  the 

distance  between  two  points.  Using  the  metric  tensor  and  (3.6),  rewrite  (3.9)  in  a 
different  coordinate  system: 

ds2  =  grj,dx'  dxj  =  gjjdx,dxJ  =  gj.A^A^dx'  dxJ  (3-8) 

From  the  above  equation,  it  can  be  seen  that  the  metric  tensor  itself 
changes  under  a  transformation: 

gir=KA],gii  (3.9) 

Typically  the  metric  tensor  is  ignored  in  Cartesian  coordinates  because  it  equals 
the  identity  tensor,  but  this  property  will  be  important  later  on  in  determining  the 
metric  tensor  of  the  transformed  space. 

Writing  (3.11)  out  in  matrix  notation  and  letting  g  =  G  and  gr.,  =  G'  it  is 
seen  that 


G'  =  A  GA 


(3.10) 


with  A  denoting  the  transformation  matrix  defined  in  (3.3).  Additionally,  it  can  be 
shown: 


g., .,  =  A*' A J.,S..  =  5, 

J  1  J  V  1  J 


(3.11) 


Also  note  that  the  inverse  of  the  metric  tensor  matrix  is  denoted  by  placing 
the  indices  in  the  superscript  position. 


«9=(G") 


(3.12) 


The  metric  tensor  not  only  characterizes  the  length  in  arbitrary 
coordinates  but  also  volume.  The  standard  Cartesian  volume  element  is 
transformed  using  the  transformation  matrix  determinant  as  indicated  below: 


dV  =  detA|^F 


(3.13) 
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Using  the  matrix  representation  (3.12)  it  can  be  seen 

g'  =  (det  A)2  g  (3.14) 

with  gand  g 'denoting  the  determinants  of  their  respective  metric  tensors.  This 
then  implies 

dV  =  Jg'dV'  (3.15) 

The  primes  can  then  be  dropped  and  the  volume  element  is  always  described  by 
■sfgdV  in  either  Cartesian  or  curved  coordinates. 

d.  Vector  Products 

For  describing  electromagnetics  utilizing  the  arbitrary 
representation  from  above,  the  vector  product  needs  to  be  defined.  Keeping  in 
mind  that  the  vector  product  is  antisymmetric,  or  UxV  =  -VxU ,  and  the  vector 
products  of  the  Cartesian  basis  vectors  are  cyclic,  define  the  Levi-Civita  tensor 

e*=[  ijk]  (3.16) 

with  [ijk]  being  the  permutation  symbol 

+1  if  i,j,k  is  an  even  pennutation  of  1,2,3 
[ijk]  =  <  -1  if  i,j,k  is  an  odd  permutation  of  1,2,3  (3.1 7) 

0  otherwise 


A  transformation  can  then  be  applied  to  (3.18)  using  the  rules  defined  above. 
This  gives  a  definition  for  the  Levi-Civita  tensor  in  arbitrary  coordinates 


=  J'k'=  AfAfAf  [Uk]  =  det(A)  )[z'y 'k']  =  1^1 


(3.18) 


which  utilizes  the  Leibniz  formula  for  the  determinant  of  A1’ .  Also,  det(A)is  the 
determinant  of  the  transformation  matrix  A), which  is  the  inverse  of  A)'.  Then 
from  (3.16)  the  determinant  is 
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det(A)  =  ±y[g'  (3.19) 

With  the  negative  sign  utilized  if  the  transformation  changes  the  handedness  of 
the  coordinate  system  and  the  positive  sign  otherwise. 

Finally,  the  definition  of  the  Levi-Civita  tensor  in  arbitrary 
coordinates  is 


dik 


(3.20) 


This  symbol  is  required  to  compute  the  vector  product  in  arbitrary  coordinates: 

U xV  =eiJk  UjVkei  (3.21) 

where  et  is  the  i-th  basis  vector  in  the  coordinate  system. 


e.  Divergence  and  Curl 

Using  the  above  notation,  the  divergence  in  three  dimensions  and 
curl  can  now  be  defined  for  an  arbitrary  coordinate  system.  While  skipping  a 
great  deal  of  important  mathematics  in  their  derivations,  it  will  suffice  to  simply 
show  what  they  are  and  move  on  to  Maxwell’s  equations. 

The  divergence  of  a  vector  field  can  be  shown  to  be 


V-V 


1  dJgT 
Jg  dx‘ 


(3.22) 


The  curl  of  a  vector  field  V  is  defined  as 

(V  x  V),  =  \}Jk]-^T  =  [Vk]djEk 

where  d i  =  d!dxj . 


(3.23) 


2.  Maxwell’s  Equations  and  Transformation  Optics 

Since  TO  deals  with  curvilinear  coordinates,  introducing  the  above 
notation  takes  care  of  the  covariance  properties  in  such  coordinate  systems. 
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Using  the  above  identities,  Maxwell’s  equations  in  empty  space  (2.1)  through 
(2.4)  can  be  rewritten  in  an  arbitrary  coordinate  system: 


1  dJgE'  _  p 

yfg  dx‘  £o 

Vg  dx‘ 


8E‘ 


8jBk  =  noso  ——  +  noj 


(3.24) 


Rewriting  (3.26)  by  placing  all  the  indices  in  the  lower  position 
through  use  of  the  metric  tensor£;.  =  g,E\  utilizing  (3.22),  and  rearranging 
slightly  yields 


d'lgg''  Ej  _  yfgP 

dx'  so 

d4ggijHj  _  Q 

dx‘ 


1  d{±4ggjHj ) 
Mo  St 


=  U  £ 

•  o  o 


e(±4gg',Ei) 


dt 


+  M 


4gf 


(3.25) 


Examining  the  above  equations  closely,  it  is  seen  that  the  form- 
invariant  Maxwell’s  vacuum  equations  resemble  the  macroscopic  Maxwell 

equations  in  dielectric  media,  keeping  in  mind  that  B  =  — 

Mo 
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it  can  be  seen  that  the  transformed  susceptibility  tensors  are  related  to  the 
transformation  itself.  In  other  words 

z!=z".=-S,±Jgg,\  (3.30) 

The  above  equations  allow  a  change  in  coordinates  to  be 
expressed  through  a  fictitious  variation  of  the  constitutive  parameters  and  source 
terms  without  needing  to  change  the  form  of  Maxwell’s  Equations  with  every 
different  coordinate  system.  This  transformation  is  implemented  through 
transformation  media. 
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a.  Transformation  Media 

Transformation  media  implement  coordinate  transformations  in 
Maxwell’s  equations.  The  way  to  think  of  this  is  to  imagine  two  different  spaces,  a 
virtual  space  and  a  physical  space.  The  virtual  space  is  a  flat  Cartesian  or 
electromagnetic  space  .  A  non-trivial  coordinate  transformation  is  then  applied  to 
the  virtual  space,  which  defines  the  transformation  medium.  This  transformation 
is  interpreted  in  the  physical  space  with  the  medium  implementing  the  coordinate 
transformation,  showing  the  ray  trajectories. 


Figure  1.  Transformation  media  implement  coordinate  transformations.  The  left 
figure  shows  the  Cartesian  grid  of  electromagnetic  space  that  is  mapped 
to  the  curved  grid  of  physical  space  through  the  coordinate  transformation 
implemented  by  A'r  in  the  right  figure.  The  physical  coordinates  enclose  a 
hole  that  is  made  invisible  in  electromagnetic  space  (where  it  shrinks  to 
the  point  indicated  there).  Consequently,  a  medium  that  performs  this 
transformation  acts  as  an  invisibility  device.  [4] 


The  transformation  is  implemented  by  performing  the  coordinate 
transformation  on  the  constitutive  relations.  Remembering  (3.30),  the 
transformed  equation  takes  the  form 

(3.30) 

As  mentioned  earlier,  since  both  the  virtual  and  physical  spaces  are  Cartesian, 
the  metric  tensor  is  the  identity  matrix  and  it’s  determinant  equals  one.  Equation 
(3.33)  is  rewritten  in  matrix  notation. 
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(3.30) 


AAr 

8  =  Ll  = - 

det  A 

The  above  equation  nicely  deals  with  the  +/-  in  the  tensor  equation  with  the 
determinant.  Now,  using  (3.32)  as  the  guide,  the  first  order  susceptibility  tensor  is 
calculated. 


X  =  s-l  = 


det  A 


(3.30) 


b.  Impedance  Matching 

The  assumption  in  the  analysis  is  that  the  original  medium  is  free 
space.  Since  it  is  important  that  the  metamaterial  be  impedance  matched  to  free 
space  to  eliminate  reflections,  the  relative  permittivity  and  permeability  must  be 
equal.  The  impedance  of  free  space  is 


%  = 


=  120/r~377Q 


(3.30) 


Impedance  in  a  dielectric  is  defined  as 


lMi+zJ  _120,  l(l+zJ 
U  y  s0(l  +  ze)  i(l  +  Ze) 


(3.30) 


Since  materials  with  matched  impedance  are  desired,  this  implies  that  the  radical 
on  the  right  side  of  (3.37)  equals  1  or 

(1  +  Ze)  =  (l  +  Zj  (3-30) 


The  above  equation  defines  the  condition  that  the  magnetic 
susceptibility  equals  the  electric  susceptibility,  at  least  to  first  order.  Since  the 
second  order,  nonlinear  electric  susceptibility  is  found  through  Miller’s  Rule,  as 
an  ansatz  the  second  order,  nonlinear  magnetic  susceptibility  will  be  taken  to  be 
equal  to  its  electric  counterpart  as  well. 


26 


C.  NONLINEAR  FORMULATION 

Based  on  the  above  arguments,  a  formulation  for  transformed  first  order 
susceptibility  is  now  found.  From  this  approximation,  a  second  order 
susceptibility  using  Miller’s  Rule  is  derived.  With  this  second  order  susceptibility, 
the  nonlinear  polarization  field  and  magnetization  of  the  dielectric  metamaterial  is 
found  and  are  used  as  source  terms  for  a  standard  electrostatic  problem  with  the 
appropriate  susceptibility. 

F?m(r)  ~  Z[%tE/(r)El!(r) 

Adding  the  polarization  field  and  the  applied  electric  field  together,  the 
displacement  field  is  found  and  analyzed  to  determine  the  radial  component  at 
the  inner  boundary  of  the  medium.  A  radial  component  implies  that  there  are 
fields  within  the  cavity.  The  field  caused  by  these  nonlinear  response  terms  is  the 
“error  field.” 
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IV.  CYLINDRICAL  CLOAK  APPLICATION 


A  cylindrical  cloak  is  modeled  below  to  demonstrate  a  potential  CDEW 
shield  design.  The  cylindrical  cloak  is  chosen  because  it  is  relatively 
straightforward  mathematically  while  still  demonstrating  the  concepts  effectively. 
This  design  applies  a  coordinate  transformation  where  the  new  coordinate 
system  has  had  the  origin  stretched  out  to  a  radius  a  =  Ira  ,  with  the  outer  radius 
b  =  \.5m  providing  finite  spatial  dimension.  Only  the  electrical  field  is  examined. 
The  magnetic  field  would  exhibit  identical  behavior  due  to  the  impedance 
matching  requirement. 


(A) 

/N 

/ 

N 

V, 

'  > 

/ 

> 

Figure  2.  Cylindrical  cloak  example.  The  thick  blue  line  shows  the  path  of  the  same 
ray  in  (A)  the  original  Cartesian  space,  and  under  two  different 
interpretations  of  the  electromagnetic  equations,  (B)  the  topological 
interpretation  and  (C)  the  materials  interpretation.  The  position  vector  x  is 
shown  in  both  the  original  and  transformed  spaces,  and  the  length  of  the 
vector  where  the  transformed  components  are  interpreted  as  Cartesian 
components  is  shown  in  (C).  [11] 


The  MATLAB  script  file  included  in  the  appendix  calculates  the  second 
order  nonlinear  susceptibility  for  a  point  contained  within  the  transformation 
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medium.  Since  by  design,  the  transformation  and  susceptibility  are  position 
dependent,  there  exists  a  tensor  field  with  a  different  tensor  for  every  point  within 
the  medium.  It  would  be  tedious  to  list  every  result  outputted  by  the  script  file. 
Instead,  the  calculations  are  outlined  in  broad  terms  and  the  results  are  then 
shown  and  analyzed. 

A.  CALCULATIONS 

The  tools  outlined  above  now  allow  the  calculation  of  the  second  order 
susceptibility  directly  from  the  first  order,  transformed  susceptibility.  Below  are 
listed  the  steps,  in  order,  to  perform  the  calculation. 

1.  Transformation  Matrix 

The  first  step  is  to  determine  the  transformation  matrix  needed  to  move 
from  the  primed,  standard  Cartesian  coordinate  system  in  electromagnetic  or 
virtual  space,  to  the  unprimed  system  in  physical  space  which  contains  the 
medium.  In  the  unprimed  system,  the  origin  is  expanded  to  a  radius  a.  An  outer 
radius  b  gives  the  transformation  finite  spatial  extent.  The  designed  material 
transforms  the  radial  component  yet  preserves  the  angle  and  vertical 
components  according  to 

b-a  , 

r  = - r  +a 

b 

</>  =  </>'  (3.30) 

z  =  z' 

When  r'  =  0  in  the  primed  virtual  coordinate  system,  the  radius  is  “a”  in  the 
unprimed  system,  as  expected.  Additionally,  when  r'  =  b,  the  radius  now  equals 
“b”  in  the  unprimed  system. 

Utilizing  the  standard  equations  relating  Cartesian  to  cylindrical 
coordinates,  it  is  found  that 
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By  expressing  both  the  primed  and  unprimed  system  in  Cartesian  systems,  the 
metric  tensor  goes  to  the  identity  tensor  and  the  calculations  are  simplified. 

Using  the  standard  formulation  for  the  transformation  matrix  (3.3),  the 
Cartesian  transformation  matrix  from  the  unprimed  to  the  primed  system  is  [13] 


The  above  is  written  in  cylindrical  coordinate  notation  to  save  space  with 

_  dr  _b-a 


2.  First  Order  Response 

The  linear  displacement  field  is  found  using 

DIJ  =  Eij  (3.32) 

utilizing  (3.34)  for  the  permittivity.  The  permittivity  of  free  space  is  included  to 
provide  units.  The  first  order  susceptibility,  necessary  for  finding  the  second  order 
susceptibility,  is  found  using  (3.35). 

3.  Second  Order  Response 

Plugging  the  first  order  susceptibility  into  Miller’s  Rule  (2.48)  and 
performing  the  sum  will  yield  a  second  order  susceptibility  for  any  given  point 
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within  the  confines  of  the  medium,  i.e.  r  <  b  .  This  is  a  rank  3  tensor  and  thus  has 
27  elements. 

For  the  nonlinear  case,  the  nonlinear  displacement  field  is  then  found  by 
using  the  second  order  susceptibility  and  applying  it  to  (2.21),  obtaining  the 
nonlinear  permittivity.  This  permittivity  is  then  inserted  into  the  equation  for  the 
displacement  field  (2.6)  to  obtain  the  nonlinear  response. 

B.  RESULTS 

In  order  for  the  cavity  to  be  free  of  electric  fields,  the  normal  component  of 
the  displacement  field  must  equal  zero  at  the  inner  surface  of  the  dielectric 
media.  Strictly  speaking,  the  difference  between  the  normal  components  on  the 
inner  dielectric  media  boundary  equals  the  surface  charge  density  sigma. 

D1±  -  D2±  =  a  (3.32) 

Since  it  is  a  dielectric  and  not  a  conductor,  the  surface  charge  density  is  taken  to 
be  zero.  This  defines  the  boundary  conditions  necessary  to  shield  the  cavity 
against  incoming  radiation  fields.  Another  way  of  expressing  this  is  that  the  radial 
component  of  D  at  the  inner  surface  falls  to  zero,  or  that  the  dot  product  of  the 
vector  with  its  position  vector  equals  zero,  implying  they  are  orthogonal. 

For  the  below  examples,  the  applied  electric  field  equals  100,000  V/m  or 
lOOkV/m,  which  is  a  large  electric  field  strength.  While  the  intensity  of  the  applied 
field  certainly  affects  the  field  strengths  in  the  solution,  it  does  not  affect  the 
qualitative  behavior  of  the  results.  The  area  of  interest  is  in  the  nonlinear  realm, 
but  the  linear  response  is  included  for  completeness  and  to  verify  that  the 
method  of  calculation  is  correct. 

1.  Linear  Response 

The  linear  response  obtained  from  the  MATLAB  script  falls  in  line  with 
what  is  to  be  expected  from  a  standard  transformation  optics  design  technique. 
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Below  is  a  plot  of  the  coordinate  system  and  how  it  deforms  under  the 
transformation. 


Linear  D-Field  (C/m2),  Applied  E  Field  in  +y  Direction 


Figure  3.  Linear  response  to  an  applied  electric  field  oriented  in  the  +y  direction. 

The  tangential  behavior  of  the  displacement  field,  shown  in  blue,  implies 
shielding  of  the  cavity.  The  streamlines  are  added  to  illustrate  the 
geometry  of  the  solution. 

Each  streamline  is  a  line  of  constant  x  and  demonstrates  how  the  medium 
transforms  space.  As  mentioned  previously,  the  electromagnetic  fields  are  taken 
to  be  attached  to  the  coordinate  system,  and  therefore  it  is  expected  that  the 
displacement  field  vectors  follow  these  lines. 

Figure  4  shows  the  displacement  vector  field  coincident  with  the 
transformed  coordinate  lines  from  the  previous  figure.  The  displacement  field 
exhibits  tangential  behavior  with  the  inner  surface  of  the  medium,  implying  zero 
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radial  components  and  thus  shielding  of  the  linear  fields  within  the  cavity.  A  close 
up  view  in  Figure  5  shows  that  the  displacement  field  vector  near  the  boundary  is 
indeed  tangential. 


Linear  D-Field,  Applied  E  Field  in  +y  Direction 


Figure  4.  Displacement  vector  field  transformed  by  the  linear  response  of  the 
medium.  The  tangential  behavior  of  the  D  field  is  apparent. 


34 


Linear  D-Field,  Applied  E  Field  in  +y  Direction 


Figure  5.  Close  up  of  the  linear  displacement  field  as  it  approaches  the  boundary  of 
the  medium.  As  can  be  easily  verified,  the  dot  product  of  the  field  vector 
with  its  position  vector  is  zero.  Which  means  it  is  tangent  to  the  surface  of 
the  medium.  The  [X,Y]  values  are  the  position  components,  the  [U,V] 
values  are  the  field  components.  X*U  +  Y*V  =  0. 


The  radial  component  of  the  field  line  is  zero  as  it  approaches  the  inner  medium 
boundary.  The  above  linear  response  is  expected  and  validates  the  methods 
used  to  calculate  the  fields. 

Of  note  is  the  field  strength  as  it  approaches  the  inner  boundary.  Since  the 
inner  boundary  is  the  origin  in  EM  space  blown  out  to  finite  spatial  dimension  in 
physical  space,  it  is  in  effect  a  singularity.  As  the  fields  approach  this  inner 
boundary,  the  transformed  field  strength  grows  large.  This  is  not  unexpected 
however,  as  the  first  and  second  susceptibilities  have  inverse  radial  dependence. 
As  this  radius  approaches  “a”  in  physical  space,  or  zero  in  primed  virtual  space, 
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the  fields  increase.  In  the  linear  case,  the  radial  components  remain  zero 
however  and  the  cavity  remains  shielded. 


Figure  6.  SemiLog  plot  of  the  field  strength  as  a  function  of  radius.  The  exponential 
growth  of  the  field  strength  as  it  approaches  the  inner  boundary  is  shown. 

While  the  fields  grow  large,  the  energy  remains  finite,  and  thus  does  not 
violate  conservation  of  energy.  Consider  the  energy  contained  within  a  volume, 
which  includes  the  origin,  in  virtual  space.  The  energy  density  is  mapped  to 
physical  space  by  the  transformation  medium  through  the  coordinate 
transformation  process.  The  coordinate  transformation  process  however  does 
not  induce  an  energy  flux  across  the  boundary.  Therefore,  while  the  fields  grow 
large  in  physical  space,  the  energy  remains  constant  and  finite. 
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2.  Nonlinear  Response 

The  nonlinear  response  of  the  medium  is  now  calculated  from  Miller’s  rule. 
In  order  to  calculate  the  nonlinear  displacement  field,  (2.6)  is  used  with  the 
nonlinear  polarization  field  found  from  (2.21).  Miller’s  coefficient  is  set  at  the 
experimentally  derived  average  value  of  4.52x1 0-2  m2C'1.  [10] 

a.  Miller’s  Coefficient  =  4.52x1 0-2  m2C'1 

Below  is  the  result  of  the  calculations  with  Miller’s  coefficient  equal 
to  the  accepted  average,  experimentally  derived  value.  [10]  At  first,  it  is  difficult  to 
ascertain  the  radial  component  of  the  entire  medium  since  the  symmetry  is  odd. 
The  symmetry  differs  from  that  of  the  linear  case,  but  can  be  understood  when 
one  considers  that  the  second  order  susceptibility  is  a  cubic  function  rather  than 
a  direct  product.  The  vectors  in  the  second  and  fourth  quadrants  that  are  scaled 
to  a  point  which  can  be  seen  appear  to  be  tangential  and  in  the  opposite  direction 
to  the  linear  case. 
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Nonlinear  D-Field,  Applied  E  Field  in  +y  Direction 


Figure  7.  Displacement  field  with  Miller’s  coefficient  equal  to  4.52x1 0-2  m2C'1.  There 
appears  to  be  a  slight  radial  component  as  the  medium  inner  boundary  is 
approached  in  the  first  quadrant.  The  box  is  the  area  displayed  in  the 
following  figure. 


Zooming  in  on  the  area  boxed  in  the  figure  demonstrates  that  there  are  radial 
fields  present.  Below  is  a  close  up  view  of  a  vector  in  the  1st  quadrant  closest  to 
the  boundary  with  phi  equal  to  45  degrees. 
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Nonlinear  D-Field,  Applied  E  Field  in  +y  Direction 


Figure  8.  Close  up  of  media  boundary  with  Miller’s  coefficient  =  4.52x1 0-2  rr^C'1. 

The  radial  component  of  the  displacement  field  is  apparent.  Taking  the  dot 
product  of  the  field  with  its  position  vector  yields  a  radial  component  of 
84.95  C/m2. 


It  appears  that  as  the  inner  boundary  is  approached  in  the  first  and  third 
quadrant,  the  displacement  field  magnitude  increases  and  the  field  lines  rotate  to 
a  non-tangential  orientation.  The  strength  of  the  radial  component  in  the  region 
increases  as  the  inner  medium  boundary  is  approached,  but  using  the  same 
energy  density  arguments  as  in  the  linear  case,  the  energy  remains  finite.  It  is 
also  apparent  from  the  below  figure  that  the  radial  component  varies  as  a 
function  of  phi  for  a  constant  radius  as  a  result  of  the  symmetry. 
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SemiLog  Radial  D-Field  Strength 


Figure  9.  D  field  strength  as  a  function  of  radius  for  several  different  values  of  phi. 

The  radial  component  of  the  field  increases  by  orders  of  magnitude  in 
each  case,  demonstrating  an  inability  of  the  medium  to  shield  the  cavity 
from  nonlinear  effects. 


The  two  plots  in  Figure  9  demonstrate  that  the  radial  component  of  the 
displacement  field  at  the  inner  medium  boundary  is  not  zero,  implying  a  field 
within  the  cavity.  This  radial  component  is  due  to  the  polarization  field  blowing  up 
near  the  singularity  and  dominating  the  vector  sum  in  (2.6),  or 

D  =  soE  +  P  (3.33) 

A  closer  examination  of  the  P  field  in  this  case  illustrates  the  behavior.  As  the 
inner  boundary  is  approached,  the  polarization  field  grows  at  a  much  greater  rate 
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than  the  applied  electric  field.  Thus,  the  polarization  field  dominates  in  the  above 
sum  and  the  radial  component  grows  larger.  This  implies  that  the  inner  cavity  is 
not  shielded. 


x  1  q15  SemiLog  y  of  Field  Strength 


Radius  (m) 

Figure  10.  A  semilog  plot  of  the  polarization  and  electric  fields.  The  polarization  field 
grows  at  a  much  greater  rate  as  it  approaches  the  singularity,  causing  the 
D  field  calculation  to  be  dominated  by  the  P  field. 
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V.  CONCLUSIONS 


A.  FINDINGS 

The  above  analysis  demonstrates  a  radial  component  of  the  displacement 
vector  near  the  inner  boundary  of  the  medium  when  dealing  with  the  nonlinear 
response  as  described  by  Miller’s  rule.  This  implies  the  presence  of  a  field  within 
the  cavity.  Whether  this  field  is  strong  enough  to  cause  damage  will  of  course  be 
dependent  on  what  is  contained  within  the  shield  and  its  radiation  tolerance,  as 
well  as  the  intensity  of  the  incoming  fields. 

If  Miller’s  rule  is  taken  to  be  correct,  then  the  process  of  determining  the 
first  order  susceptibility  through  TO  in  the  cylindrical  cloaking  case  necessarily 
implies  an  inability  to  shield  the  cavity  from  nonlinear  effects.  In  other  words,  the 
process  of  using  Miller’s  rule  to  derive  the  susceptibility  from  a  first  order 
susceptibility  found  through  a  coordinate  transformation  destroys  that 
transformation  and  therefore  the  geometry  of  the  cloak  for  the  second  order 
response. 

It  is  also  possible  that  Miller’s  rule  cannot  derive  the  correct  second  order 
susceptibility  when  dealing  with  TO.  As  mentioned  earlier,  Miller’s  rule  gives 
qualitative,  not  quantitative  results.  It  is  possible  that  Miller’s  rule  is  simply 
inapplicable  when  dealing  with  TO  and  that  the  second  order  susceptibility  is  not 
accurately  described  by  it.  If  this  is  the  case,  then  an  alternate  method  needs  to 
be  derived  in  order  to  define  the  second  order  susceptibility.  Perhaps  defining  the 
second  order  susceptibility  through  TO  and  then  working  backwards  through 
Miller’s  rule  to  obtain  the  first  order  susceptibility  would  be  a  good  technique.  This 
is  left  for  future  work. 

1 .  Inner  Boundary  as  a  Singularity 

The  inner  boundary  becomes  a  singularity  in  both  the  linear  and  nonlinear 

cases.  In  the  linear  case,  while  the  field  still  blows  up,  it  remains  tangential.  This 

is  not  true  for  the  nonlinear  case.  In  the  nonlinear  case,  the  polarization  field 
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blows  up  as  it  approaches  the  boundary,  causing  the  displacement  field  to  exhibit 
a  radial  component  in  certain  regions.  The  nonlinear  second  order  susceptibility 
causes  a  rotation  of  the  polarization  vector  field.  In  turn,  when  that  vector  field 
blows  up  approaching  the  medium  the  radial  component  of  the  displacement  field 
blows  up.  A  boundary  condition  setting  the  field  to  zero  at  the  inner  medium  is 
not  applicable  because  the  field  at  that  point  is  desired.  As  mentioned  previously 
however,  the  energy  density  remains  finite. 

This  condition  has  a  similar  result  documented  in  the  literature  when  the 
problem  is  taken  out  of  the  quasi-static  limit.  [13]  For  perfect  invisibility,  rays 
traversing  the  medium  need  to  traverse  the  transformed  physical  space  in  the 
same  amount  of  time  as  they  traverse  the  flat  Cartesian  electromagnetic  virtual 
space,  otherwise  distortions  would  be  introduced.  This  means  that  rays 
straddling  the  inner  boundary  need  to  traverse  this  distance  in  the  same  amount 
of  time  as  it  takes  for  them  to  traverse  the  point  in  electromagnetic  virtual  space. 
But  the  point  in  electromagnetic  space  has  zero  spatial  extent,  meaning  that  the 
rays  need  to  be  travelling  infinitely  fast.  The  conclusion  drawn  from  this  in  the 
literature  is  that  complete  invisibility  is  not  possible.  [4]  Complete  invisibility  is  not 
necessary  in  this  application  however,  and  further  study  would  need  to  be  done 
to  determine  if  the  fact  that  the  rays  are  unable  to  traverse  infinitely  fast  implies 
that  the  rays  are  not  completely  shielded  from  the  inner  cavity,  in  both  the  linear 
and  nonlinear  cases. 

B.  SPACE  APPLICATION  FEASIBILITY 

Metamaterials  themselves  have  no  intrinsic  reason  why  they  cannot  be 
deployed  on  spacecraft.  The  materials  thus  far  experimented  with  have  been 
constructed  of  metals  and  plastics,  but  the  designed  electrical  properties  of  the 
materials  are  what  is  important  and  not  the  actual  materials.  Ceramics  are  also 
being  investigated.  In  fact,  research  is  currently  underway  using  metamaterials 
and  TO  to  design  antennas  and  electrical  routing  equipment  that  would  require 
far  less  volume  and  power  to  operate  in  space.  [14]  [15] 
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The  requirement  of  complete  cloaking  is  also  not  of  concern  for  space- 
based  applications.  A  semicircular  or  flat  plate  placed  in  the  expected  direction  of 
a  directed  energy  weapon  design  to  divert  the  fields  away  would  prove  to  be 
sufficient.  In  addition,  shielding  of  only  critical  components  as  a  design 
philosophy  is  also  feasible. 

In  general,  there  are  no  outright  reasons  why  metamaterials  and  TO 
design  techniques  would  not  be  suitable  for  spacecraft  deployment.  There  are  no 
power  requirements,  weight  is  determined  by  the  materials  themselves  but  is  not 
limited  to  massive  or  fragile  materials,  and  the  spatial  extent  or  size  is  not 
necessarily  a  concern  depending  on  the  design. 

C.  CONCLUSIONS  AND  FUTURE  WORK 

Whether  or  not  TO  offers  a  way  to  shield  against  nonlinear 
electromagnetic  fields  remains  to  be  seen,  but  this  work  is  a  first  step  in  that 
direction.  Taking  the  problem  out  of  the  quasi-static  arena  and  experimentation 
with  high  strength  fields  is  necessary  to  increase  understanding.  This  is  left  as 
future  work. 

Additional  future  work  would  include  a  more  rigorous  derivation  of  Miller’s 
coefficient.  While  the  coefficient  has  been  experimentally  determined,  very  little 
progress  has  been  made  in  deriving  a  theoretical  framework  for  why  this  is  so. 
[10]  A  quantum  mechanical  or  statistical  mechanical  approach  may  help  to  shed 
light  of  this  problem. 

Finally,  as  mentioned  previously  however,  perfect  cloaking  is  not  the 
desired  end  state  of  this  research,  rather  the  elimination  or  diminution  of  the 
fields  internal  to  the  cavity  to  shield  it  from  damage.  The  presence  of  normal  field 
components  demonstrated  above  implies  that  the  nonlinear  fields  will  “leak”  into 
the  cavity.  The  strength  of  these  fields  is  what  next  needs  to  be  determined.  The 
calculations  above  say  it  would  be  an  infinite  field,  but  this  is  obviously 
nonphysical.  Future  work  needs  to  be  conducted  to  determine  the  strength  of 
these  fields. 
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Transformation  optics  and  the  use  of  metamaterials  as  transformation 
media  is  a  nascent  yet  exciting  development  in  electrodynamics.  While  the  field 
of  linear  transformation  optics  has  been  developed  to  a  point  where  real 
materials  are  producing  results  in  the  lab,  the  nonlinear  effects  are  still  a  wide- 
open  subject  for  theoretical  as  well  as  experimental  studies. 
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APPENDIX.  MATLAB  CODE  USED  FOR  CALCULATIONS 


%{ 

THESIS  SCRIPT 

LCDR  MATTHEW  DEMARTINO 

14DEC2012,  NAVAL  POSTGRADUATE  SCHOOL 

NOTES 

THIS  SCRIPT  MAKES  SEVERAL  FUNCTION  CALLS  AND  CALCULATIONS  NECESSARY  TO 
PLOT  THE  LINEAR  AND  NONLINEAR  DISPLACEMENT  FIELDS  FOR  A  CYCLINDRICAL 
CLOAK  OF  OUTER  RADIUS  1.5  M  AND  INNER  RADIUS  1M.  THE  APPLIED  ELECTRIC 
FIELD  IS  lOOkV/m  ORIENTED  IN  THE  POSITIVE  Y  DIRECTION 

%} 

clear  all 
close  all 
clc 

I  =  [1  0  0;0  1  0;0  0  1];  % Identity  matrix 

%Define  the  virtual  prime  space  extent  and  discrete  step 
step  =  .01; 
minstart  =  -1.5; 
maxstart  =  1.5; 

rng  =  (maxstart-minstart) /step; 

[X,Y]  =  meshgrid (minstart : step :maxstart) ;  %X'  and  Y'  in  flat  EM  space 

%Define  the  transformation  medium 
a  =  1;  %inner  radius  (m) 

b  =  1.5;  %outer  radius  (m) 

z  =  0 ; 

eo  =  8.854e-12;  %permittivity  of  free  space,  SI  units 

%Initialize  electric  field  vectors  for  each  point  x,y  --  P  =  X  Eo(jhat) 
EX  =  zeros (size (X) ) ; 

EY  =  zeros (size (Y) ) +le6;  %E  =  0 [x]  +  100000 [y] 

%Initialize  other  data  structures 

%P  fields,  "2"  =>  2nd  order,  "N"  =>  normalized 

PX  =  zeros (size (X) ) ; 

PY  =  zeros (size (Y) ) ; 

PX2  =  zeros (size  (X)  )  ; 

PY2  =  zeros (size (Y) ) ; 

PX2N  =  PX2 ; 

PY2N  =  PY2 ; 

Z  =  zeros (size (X) ) ;  %Used  for  contour  plots  if  desired 

%D  fields,  "2"  =>  2nd  order,  "N"  =>  normalized 
DX  =  zeros (size (X) ) ; 
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DY  =  zeros (size (X) ) ; 

DXN  =  zeros (size (X) ) ; 

DYN  =  zeros (size (X) ) ; 

DX2  =  zeros (size  (X)  )  ; 

DY2  =  zeros (size (X) ) ; 

DX2N  =  zeros (size (X) ) ; 

DY2N  =  zeros (size (X) ) ; 

%loops  to  step  through  each  x  and  y  coordinate  in  virtual  space  and 
transform  it  and  its  fields  to  physical  space 
for ( i  =  1 : rng) 

for(j  =  l:rng) 

rprime  =  sqrt (X ( i , j ) A2+Y ( i , j ) A2 ) ; 
r  =  a  +  rprime* (b-a) /b; 

R  =(b-a)/b;  % (dr/dr' ) 

phi  =  atan2 (Y (i, j ) , X (i, j ) ) ; 

if (r  <  b) 

%Calculate  transformation  matrix  for  point  [X ( i , j ) , Y ( i , j ) ] 
A(l,l)  =  a/ (rprime)  -  (a*X (i, j ) A2) / (rprimeA3)  +  (b-a) /b; 

A  ( 1 , 2 )  =  (-a*X (i, j ) *Y (i, j ) ) / (rprimeA3) ; 

A ( 1 , 3 )  =  0; 

A  (2 , 1 )  =  (-a*X  (i,  j  )  *Y  (i,  j  )  )  /  (rprimeA3)  ; 

A(2,2)  =  a/ (rprime)  -  (a*Y (i, j ) A2) / (rprimeA3)  +  (b-a) /b; 

A (2 , 3 )  =  0; 

A ( 3 , 1 )  =  0; 

A ( 3 , 2 )  =  0; 

A ( 3 , 3 )  =  1; 

%calculate  the  first  order  susceptibility 

eps  =  eo* (A*A' ) . / ( (r/rprime) *R) ; %f irst  order  trans  epsilon 
Xi  =  eps/eo-I;  %first  order  susceptibility 

^transform  the  electrical  fields  to  their  physical  space 
% (unprimed)  values 

EX ( i ,  j  )  =  (EX ( i , j ) *A (1,1) +EY (i , j ) *A (1,2) ) ; 

EY  ( i , j )  =  (EX (i,j)*A(2,l) +EY (i , j ) *A (2,2) ) ; 

%explicit  matrix  vector  product  for  first  order  polariation 
PX(i,j)  =  eo*  (Xi ( 1 , 1 )  *  EX ( i , j )  +  Xi(l,2)  *  EY(i,j)); 

PY(i,j)  =  eo* (Xi (2,1)  *  EX ( i , j )  +  Xi(2,2)  *  EY(i,j)); 


%second  order  susceptibility 
Xi2  =  SecOrdSuscept2 (Xi) ; 

%2nd  order  polarization  explicit  calculation (2D) 
PX2 (i, j )  =  eo* ( (Xi2 (1,1,1) *EX(i,j) *EX(i,j) )  + 

(Xi2 (1,1,2) *EX(i, j) *EY (i, j) )  +  (Xi2 ( 1 , 2 , 1 ) *EY ( i , j ) *EX ( i , j ) )  + 
(Xi2 (1,2,2) *EY(i, j) *EY(i, j) ) ) ; 

PY2 (i, j )  =  eo* ( (Xi2 (2,1,1) *EX(i,j) *EX(i,j) )  + 

(Xi2 (2,1,2) *EX(i,j) *EY(i,j) )  +  (Xi2 (2,2,1) *EY (i, j ) *EX (i, j ) )  + 
(Xi2 (2,2,2) *EY(i, j) *EY(i,  j) ) )  ; 

PX2N ( i , j )  =  PX2 (i, j) /sqrt(PX2 (i, j) A2+PY2 (i, j) A2) ; 
PY2N ( i , j )  =  PY2 (i, j) /sqrt(PX2 (i, j) A2+PY2 (i, j) A2) ; 
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%calculate  D  fields 

%"N"  =  normalized,  "2"  =  2nd  order 

DX  ( i , j )  =  eps (1, 1) *EX (i, j ) +  eps (1,2) *EY  ( i , j ) ; 

DY ( i , j )  =  eps (2, 1) *EX  (i, j ) +  eps (2, 2) *EY  (i, j )  ; 

DXN  ( i , j )  =  DX(i,  j) /sqrt (DX(i, j) A2+DY(i, j) A2)  ; 

DYN ( i , j )  =  DY(i, j) /sqrt(DX(i, j) A2+DY(i, j) A2) ; 

DX2  ( i , j )  =  (eo*EX (i, j ) +PX2  (i, j ) )  ; 

DY2  ( i , j )  =  (eo*EY (i, j ) +PY2  (i, j ) ) ; 

DX2N ( i , j )  =  DX2  (i,  j  ) /sqrt  (DX2  (i,  j  )  A2+DY2  (i,  j  )  A2)  ; 
DY2N ( i , j )  =  DY2 (i, j ) /sqrt (DX2 (i, j ) A2+DY2 (i, j ) A2) ; 

%transform  the  points 
X(i,j)  =  r*cos(phi); 

Y  (i, j )  =  r*sin (phi) ; 


end 

end 

end 


%Plotting  functions 
%PlotLinear ; 

%PlotNonLinear ; 

function  [X2]  =  SecOrdSuscept (Xi ) 

%This  function  takes  as  input  a  first  order  susceptibility 

%and  returns  a  second  order  susceptibility  based  upon  Miller' s  Rule 

%Calculate  the  second  order  susceptibity 

X2  =  zeros ( 3 , 3 , 3 ) ;  %rank  3  tensor  =>  27  elements 

for ( i  =  1:3); 

for ( j  =  1:3) 

for ( k  =  1:3) 

%Perform  Miller' s  rule  sum 

for (p  =  1:3) 

for (q  =  1:3) 

for (r  =  1:3) 

X2(i,j,k)  =  X2 ( i , j , k)  +  Xi (i,p) *Xi (q, j ) *Xi (r, k) 

end 

end 

end 

end 

end 

end 

K  =  4.52e-2;  %Miller' s  coefficient  in  mA2/C 
X2  =  -K*X2 ; 
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